function [v_partial_dt, v]  = time_derivative_v(v_fitted,  model_sol,   par,    hat_kappa0 )
% this function returns the matrix of v_tilde derivatives, after updating
% from the last round v_tilde functions.   v_fitted  is a scatter interpolant object. 
w_grid = par.w_grid;  lambda_grid=par.lambda_grid;
tau=par.tau;   pi=par.pi;  
if( isfield(par, 'kappap_multiplier') )
    kappap_multiplier = par.kappap_multiplier;
else
    kappap_multiplier = 1;
end

if( par.model_name~="benchmark" )   % for non-benchmark models. 
    w = model_sol.w_matrix; 
    lambda=model_sol.lambda_matrix; 
    lambda_tau = 1/tau - pi*lambda;
    % get the derivatives
    v = v_fitted( w,  lambda );
    v_partial_w = v;  v_partial_lambda=v;    v_partial2_w=v;
    for( j = 1:par.N_lambda)
        v_partial_w(:,j) =   derivative( w_grid, v(:,j) ) ;
        v_partial2_w(:,j) =  derivative( w_grid, v_partial_w(:,j) ) ;
    end

    for( i = 1:par.N_w)
        v_partial_lambda(i,:) = derivative( lambda_grid,  v(i,:)' )  ;
    end

    kappaw = model_sol.kappaw_matrix;
    jump_lambda = par.kappa_lambda_fun(lambda);
    Delta_v = v_fitted(  w,  lambda ) - v_fitted(  w.*(1-kappaw), lambda+jump_lambda );

    muw = model_sol.muw_matrix;
    mu_lambda = par.mu_lambda_fun(lambda);
    xK = model_sol.xK_matrix;
    kappab = model_sol.kappab_matrix;
    sigmaw = model_sol.sigmaw_matrix;
    sigmaK = par.sigmaK;
    sigmap = model_sol.sigmap_matrix;
    if(hat_kappa0==0)   % Note: if hat_kappa0=0, then we are solving a safe bond pricing. Therefore, we set the jump part to zero as well. 
        kappap = 0;
    else
        kappap = kappap_multiplier*model_sol.kappap_matrix;  % use the multiplier over the fundamental kappap to price credit. 
    end
    rf = model_sol.rf_matrix;
    v_partial_dt = lambda.*pi./(1-kappab) .* ( v - (1-kappap-hat_kappa0) ) + lambda.*(1-pi)./(1-kappab).*Delta_v + lambda_tau.*(v-1) + xK.*(sigmaK+sigmap).^2.*v_partial_w.*w.*sigmaw + rf .* v -  (  v_partial_w.*w.*muw + 0.5*v_partial2_w.*w.^2.*sigmaw.^2 + v_partial_lambda.*mu_lambda  ) ;
else
    w = w_grid; lambda=mean(lambda_grid);   lambda_tau = 1/tau - pi*lambda;
    v = v_fitted( w,  lambda*ones(size(w)) );   % a column vector    
    v_partial_w = derivative( w_grid, v );
    v_partial2_w = derivative( w_grid, v_partial_w );

    muw = mean(model_sol.muw_matrix,2);
    xK = mean(model_sol.xK_matrix,2);
    kappab = mean(model_sol.kappab_matrix,2);
    sigmaw = mean(model_sol.sigmaw_matrix,2);
    sigmaK = par.sigmaK;
    sigmap = mean(model_sol.sigmap_matrix,2);
    kappap = kappap_multiplier*mean(model_sol.kappap_matrix,2);  % For robusntess, take the average. They should be the same across all lambda.
    rf = mean(model_sol.rf_matrix,2);
    
    v_partial_dt = lambda.*pi./(1-kappab) .* ( v - (1-kappap-hat_kappa0) ) + lambda_tau.*(v-1) + xK.*(sigmaK+sigmap).^2.*v_partial_w.*w.*sigmaw + rf .* v -  (  v_partial_w.*w.*muw + 0.5*v_partial2_w.*w.^2.*sigmaw.^2  ) ;
    % Finally, reshape everything
    v = v * ones( 1, par.N_lambda ) ;
    v_partial_dt =  v_partial_dt * ones( 1, par.N_lambda ) ;
end








